Chemical Reaction Dynamics within Anisotropic Solvents in Time-Dependent Fields 



Eli HershkovitfE and Rigoberto Hernandez 
Center for Computational Molecular Science and Technology 
School of Chemistry and Biochemistry 
Georgia Institute of Technology 
Atlanta, GA 30332-040^ 
(Dated: February 2, 2008) 

The dynamics of low-dimensional Brownian particles coupled to time-dependent driven anisotropic 
heavy particles (mesogens) in a uniform bath (solvent) have been described through the use of a 
variant of the stochastic Langevin equation. The rotational motion of the mesogens is assumed 
to follow the motion of an external driving field in the linear response limit. Reaction dynamics 
have also been probed using a two-state model for the Brownian particles. Analytical expressions 
for diffusion and reaction rates have been developed and are found to be in good agreement with 
numerical calculations. When the external field driving the mesogens is held at constant rotational 
frequency, the model for reaction dynamics predicts that the applied field frequency can be used to 
control the product composition. 



I. INTRODUCTION 

The stochastic or Brownian motion of a particle in 
a uniform solvent is generally well-understoodi^ The 
dynamics is less clear when the solvents respond in a 
non-uniform or time-dependent manner, although such 
problems are not uncommon. For example, the dy- 
namical properties of a suspension in a liquid crystal 
can be projected onto an anisotropic stochastic equa- 
tion of motion. "^'^ Other examples may include diffu- 
sion and reaction in supercritical liquids^ liquids next 
to the liquid vapor critical point-'-'- and growth in living 
polymerization.^ 

The flow properties of liquid crystals have gener- 
ally been analyzed from the perspective of macroscopic 
nematohydrodynamicsjifi Therein, liquid crystals have 
been classified according to the presence or absence of 
solvent. Pure liquid crystals containing no solvent are 
called thermotropic in part because they have exhibited 
strong temperature-dependent behavior. A suspension of 
nematogens (anisotropic molecules) within a simple sol- 
vent is known as a lyotropic liquid. The presence of ne- 
matogens leads to different transport properties within 
the solvent than would be seen in a pure simple liquid 
alone. The additional complexity is a result of the cou- 
pling between the velocity field and the average direction 
of the nematogens. As a result, the dynamics of a parti- 
cle in the liquid crystal is dissipated by a friction whose 
form is that of a tensor and not a scalariii The actual 
drag can be further complicated by the presence of topo- 
logical discontinuities in the liquid. To our knowledge, 
analytic solutions for the diffusion of Brownian particles 
in these general environments are not known. The situa- 
tion for a reactive solute is even less clear as no analytic 
formalism has been constructed. In the present work, we 
construct a formalism — that in some limits — fills in these 
gaps. 

One step toward understanding the dynamics in 
anisotropic liquids would thus be the development of a 
lyotropic model consisting of a Brownian particle in the 



presence of a time-dependent driven mesogen^ Another 
step toward this goal is the analytic and/or numerical 
solution of such. In the present work, the rigorous con- 
struction necessary for the first of these steps is not at- 
tempted. Instead, a naive phenomenological model de- 
scribing the dynamics in lyotropic liquids has been con- 
structed. It serves as a benchmark for the development of 
techniques useful in analyzing the dynamics of Brownian 
particles dissipated by an anisotropic solvent through a 
time-dependent friction. In particular, the lyotropic liq- 
uid is assumed to be nematic, i.e., the (calamitic) meso- 
gens are assumed to be rod-like as is the case with mineral 
moieties in water—. The mesogens are further assumed 
to be one-dimensional and rigid, and a series of additional 
simplifying assumptions have been invoked. A physical 
system rigorously satisfying all these assumptions may 
not exist, but the benchmark may still exhibit some of 
the important dynamics that has been seen in real liquid 
crystals in the presence of magnetic fields with time and 
space instabilities^ Another step toward understanding 
the dynamics in anisotropic liquids is the rigorous solu- 
tion of a thermotropic (nematic) model in which the di- 
lute Brownian particle diffuses or isomcrizcs in a solvent 
that consists exclusively of mesogens. It is based on the 
possible connection to a rotating nematic liquid system 
previously observed, and on the analytic understand- 
ing of the dynamics in nematic liquids in a few special 
casesiiiiSi For this thermotropic case, we don't attempt 
to develop a microscopic model of the friction and in- 
stead make assumptions based on the known properties 
of isotropic liquids. 



In general, the complicated microscopic dynamics of a 
subsystem coupled to a many-dimensional isotropic heat 
bath can be projected onto a simple reduced-dimensional 
stochastic equation of motion in terms of the variables of 
the subsystem alone. In the limit when the fluctuations 
in the isotropic bath are uncorrelated, the equation of 
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motion is the Langevin Equation (LE)^ 

q = p (la) 

P = -V'iq)--/p + m ■ (lb) 

where {q,p) are the position and momenta vectors in 
mass- weighted coordinates (i.e. mass equals one), V{q) 
is the system potential, 7 is the friction and ^ is a Gaus- 
sian random force due to the thermal bath fluctuations. 
The friction and the random force are connected via the 
fluctuation dissipation theorem, 

mnt2)) = ^s{h-t,), (2) 

where the average is taken over all realizations of the 
forces at the inverse temperature (3[= (fceT)"^]. The 
LE can represent the generic problem of the escape rates 
of a thermally activated particle from a metastable well 
when the thermal energy is much lower then the bar- 
rier height jSi The one-dimensional LE has been solved 
in the asymptotic limits of weak and strong friction by 
Kramersi^ A general solution for weak to intermediate 
friction was found by Melnikov and Meshkovii^ This 
result was subsequently extended to the entire friction 
range in the turnover theory of PoUak, Grabert and 
Hanngi;-^ The reactive rates for a multidimensional LE 
have been obtained exactly in the strong^^'^^i^'^i^^i^^i^^ 
and weak frictior^ limits and approximately in be- 
tween these limits through a multidimensional turnover 
theoryiS. The LE can also describe the dynamics of a sub- 
system under an applied external force, and has led to 
the observation of such interesting phenomena as stochas- 
tic resonancerSiSS resonant activation^SL^i and rectified 
Brownian motioni^Si^^iS^iS 

When the fluctuations in the isotropic bath do not de- 
cay quickly in space or in time, the dynamics are known 
to be described by the the Generalized Langevin Equa- 
tion (GLE)f3fii The activated rate expression for a parti- 
cle described by a GLE is also well-knowni^LSi Less un- 
derstood are the exact rates when the friction dissipates 
the subsystem differently at different times in a nonsta- 
tionary GLE-like equation.Si2iii22iiMi Nonetheless, the 
models developed in this work contain the flavor of this 
nonstationarity in that the LE is driven by an external 
periodic field through the friction rather than through a 
direct force on the system. Consequently the result of 
this study also provide new insight into the dynamics of 
systems driven out of equilibrium. 

The primary aim of the paper is the development of 
analytical and numerical techniques to obtain the dif- 
fusion and reaction rates of a subsystem dissipated by a 
time-dependent driven anisotropic solvent in various lim- 
its. A naive model for a nematic lyotropic liquid and its 
various underlying assumptions is presented in Sec. HTIas 
one paradigmatic example for the accuracy of the meth- 
ods described in this work. Another model based on 
an experimental system of the rotating nematic liquid 
is described briefly in Sec. Illll The anisotropic solvent is 



manifested in these models by way of a time-dependent 
friction that is externally driven. The diffusion of free 
Brownian particles dissipated by a time-dependent envi- 
ronment is described in Sec. lIVI The numerical methods 
for calculating reactions rates needed to extend the solu- 
tions of these models to include nontrivial potentials of 
mean force are presented in Sec. IV Al Analytical approx- 
imations for otherwise-rigorous rate formulas are derived 
and compared to the the numerical results in Sec. IV Bl 
A discussion of the validity of all of these approaches and 
possible applications concludes the paper in Sec. IVII 

II. A NAIVE LYOTROPIC MODEL WITH 
ROTATING EXTERNAL FIELDS 

A naive model describing a particle propagated in an 
anisotropic solvent is motivated in this section in the con- 
text of diffusive or reactive dynamics within a lyotropic 
solvent. The connection between the model and realiz- 
able lyotropic solvents is only a loose one. No attempt is 
made here to do a rigorous projection of the detailed com- 
plex modes of the lyotropic solvent onto the subsystem 
dynamics. The lyotropic liquid is assumed to consist of 
rod- like mesogens and a uniform isotropic liquid solvent. 
It is further assumed that there exists a single tagged 
motion characterized by an effective coordinate q that 
describes the subsystem — e.g., a probe particle or react- 
ing pair of particles — whose dynamics is of interest. This 
tagged motion is taken to be one-dimensional for simplic- 
ity. The effective mass rUg associated with the tagged 
subsystem is also assumed to be well separated from the 
smaller mass of the isotropic liquid, and the larger mass 
of the (anisotropic) rod-like mesogens. Consequently the 
tagged motion can be described as that of a Brownian 
particle at position q experiencing a dissipative environ- 
ment due to the interactions with the isotropic liquid and 
the mesogens. 

The model is further simplified by assuming that the 
mesogens of given concentration, c, do not interact with 
each other. This ideal-solute assumption is certainly real- 
ized at low enough concentrations that the mean spacing 
between mesogens is long compared to their effective in- 
teraction distance. (It would be easy to achieve such con- 
centrations even at relatively high concentrations if the 
interaction potentials are hard-core.) The ideal-solute 
mesogens will exhibit no orientational order in the ab- 
sence of external fields. 

In real nematic liquids there are interactions between 
the mesogens that result from cooperative forces. They, 
as well as boundary effects on the rods, are excluded 
within the model of this work. The orientation of all the 
rods is firmly fixed by a magnetic field (homogeneous 
director field) with inclination 6 relative to the y axis: 

Hx — Hq sin 9 (3a) 
Hy = ijocose* (3b) 
i?. = 0, (3c) 



3 




FIG. 1: A Brownian particle with a diameter 2a moves with 
velocity v inside a mixture of an isotropic liquid and calamitic 
mesogens. The mesogens have a length I of the same order of 
a, a negligible width and concentration, c. The mesogens are 
oriented by an external magnetic field H. The magnetic field 
is characterized by the angle 6 relative to v. Under these con- 
ditions, the Brownian particle collides with {-kFI? + 2Rl \ cos ^|) 
mesogens per unit time. 

This strong field assumption — all the mesogens will ori- 
ent uniformly in the direction of H — also ensures that 
there is no angular momentum transfer in collisions be- 
tween the mesogens and diffusing Brownian solutes. The 
environment is clearly anisotropic, and a Brownian parti- 
cle diffusing through it would experience different dissipa- 
tive forces depending on the direction of its motion. The 
suspended particle is assumed to have a spherical shape 
with a radius R. The particle velocity v{t) is restricted 
to the X direction. The number of collisions per unit 
time between the Brownian particle and the mesogens is 
simply {ttR^ + 2Rl \ cos6'|)wc. This result is illustrated in 
Fig. n Further assuming that each of the mesogens has a 
thermal distribution of velocities and noting the mass dif- 
ference between the mesogens and the environment, the 
friction force on the Brownian particle is proportional to 
the number of collisions. This gives a friction coefficient: 

7 = 7o(7ri?^ + 2i?/|cos6'|) , (4) 

where 70 is a proportionality constant characteristic of 
the system. The viscosity of the isotropic liquid leads to 
additional dissipation that is manifested as an additional 
constant term to the overall friction. However, in those 
cases when this isotropic friction is dominated by the 
friction of Eq. ^ its effect is small and insufficient to 
blur the anisotropy of the system. For further simplicity, 
therefore, in what follows, the isotropic friction due to the 
homogeneous solvent will be assumed to be zero without 
loss of generality as long as the actual isotropic friction 
is weak in this sense. 

The instantaneous inclination 9{t) has a large influence 
on the short-time dynamics of a particle whose motion is 
measured only along an initially chosen x axis. Without 



this restriction, a fixed 9 will not influence the dynamics 
because the particle motion will necessarily average over 
all directions. As a result, the inclination can be used as a 
control parameter when one measures only the dynamics 
along a specific direction but not when one is interested in 
the average diffusion or reaction of the chosen subsystem 
with respect to all directions. 

In cases when the magnetic field of Eq. O rotates with 
frequency then the Brownian particle experiences a 
friction, 

7(t) =7o(7ri?2-f 2ffl|costjt|) , (^5^ 

that is periodic in time. Including the dissipation of the 
rotating mesogens will not change this friction, but will 
add a finite temperature to the bath due to rotational dis- 
sipation. As long as this amount of heat is much smaller 
then the bath temperature, the friction in Eq. [Sjis well- 
defined and can be used as the friction entirely dissipating 
the Brownian particle. 

III. A NEMATIC MODEL WITH EXTERNAL 
ROTATING FIELDS 

While the naive model described above does capture 
some of the features of liquid crystal diffusion, it is 
nonetheless too simplistic. Experiments of pure nematic 
liquids under a rotating magnetic field^^*^ can serve to 
illustrate the possibility of solvent responses character- 
ized by time dependent viscosity. In these experiments, 
the homogeneous director field of a nematic liquid con- 
fined between two parallel glass plates was aligned in the 
plane of the plates by strong magnetic field. The mag- 
netic field was also rotated at constant velocity within 
this plane. For many of the experimental conditions, the 
nematic liquid retained uniform alignment but its homo- 
geneous director field followed the magnetic field with a 
constant phase lag. Finding an expression for the viscos- 
ity in a nematic liquid is far more complicated then for 
an isotropic liquid. It has five coefficients and depends 
on the orientation of the director, the velocity and the ve- 
locity gradient. This problem was only partially solved 
for some special cases. One case obtains the effective 
viscosity in a suspension of small particles in a nematic 
liquid.^* The key simplifications are that the small parti- 
cles are assumed to be not much larger than the nemato- 
gens and with spherical shape. The friction coefficient 
has the simple form, 

f,^a{A6,k+BcosH), (6) 

where the expression for constant coefficients A and B 
may be found in Ref. 0. A second case treats the limit 
in which the chosen particle in a nematic liquid has a 
large spherical shapeiii The resulting effective friction is 
composed of an isotropic term and an anisotropic term 
that depends on the angle between the director and the 
particle velocity. The anisotropic expression is a little 
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more complicated than Eq. but its leading order terms 
also involve sin 9 and cos 6'. In both of these cases, the 
nematic liquid is assumed to be firmly oriented by a 
strong external field and the friction force is taken to 
be much larger than the clastic forces in the nematic 
liquid. Thus the naive model described in the previous 
section does exhibit both a uniform constant term and 
an anisotropic oscillatory term that are in qualitative — 
though not quantitative — agreement with more detailed 
models. 



IV. FREE BROWNIAN DIFFUSION IN AN 
ANISOTROPIC SOLVENT 

The motion of a free Brownian particle in the time- 
dependent friction field of Eq. [S] can be described by the 
Langevin equation, 

p--7o0Wp + ^WeW , (7) 
where the time-dependent coefficients, 

Ht) = (8a) 
ipit) = a + cos{ujt) , (8b) 
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FIG. 2: The mean square displacement of a free Brownian 
particle in the naive lyotropic bath model has been obtained 
by direct integration and through the use of the analytical 
expression in Eq. 1131 at various frequencies uj of the driving 
rotating magnetic field. In the former integration method, 
100,000 trajectories were suflicient to obtain convergence. In 
the latter, the average is taken over an ensemble of 100,000 
particles starting at time t — with inclination perpendicular 
to the velocity, and overlays the results of the former within 
the resolution of the figure. 



have been chosen to describe the periodic behavior of the 
naive lyotropic model and the hydrodynamical friction 
terms in pure nematics as simply as possible. The noise 
is related to the friction by the fluctuation dissipation 
relation, 



(9) 



The strength a of the isotropic term has been chosen to 
be 1.05 throughout the illustrations in this work to em- 
phasize the anisotropic effects, but different physically- 
realizable strengths do not lead to different conclusions. 
The solution to the equation of motion [T] is 

p{t) = poexp[-7oG(t)] 

+ / dtiVXti)e(ii)e-^°^^'*^-^(*^)> , (10) 



where po and to satisfy the initial condition, po = p{to), 
and the integrated friction G{t) is defined as 



G{t) = / dti^pity 

J to 



(11) 



The velocity correlation function is readily calculated to 
be 



(pitMh)) = iexp[-7o{G(ii) + G(i2) 
-2G(min(ti,i2))}] • 



The square mean displacement of the free particle after 
time t is the double integral. 



t nt 



dtldh {p{tl)p{h)) 



to J to 



(13a) 
. (13b) 



(12) 



'to J to 

+ f dt, rdi2e[''°{«(*^)-^(*^)>] 



A similar expression was developed by Drozdov and 
Tucker'^ for the case of fiuctuations in the local density 
of supercritical solvent. The result in Eq. ^| leads to 
the diffusivity of the particle. As will be shown below, 
the diffusivity in the time-dependent environment devi- 
ates from the linear correlation known to result in the 
constant friction environment. 

In Fig. 121 the mean square displacement of a Brow- 
nian particle whose motion is measured only along an 
arbitrary one-dimensional axis is plotted as a function 
of time at various applied frequencies lo. The average 
behavior of the mean square displacement is linear with 
time, as in the constant friction regime, but it also con- 
tains fluctuations (in time) around the average whose fre- 
quency depends on the external field. It is important to 
note that the overall slope of the mean square displace- 
ment depends on the frequency w; that is, the diffusivity 
shows strong dependence on lo. Hence by changing the 
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frequency of the external field, it becomes possible to 
control the diffusivity of the Brownian particles. 

The analytical result of Eq. E| was used to check the 
accuracy of the numerical integrator employed in prop- 
agating particles in a time-periodic white noise bath. A 
fourth-order integrator was developed based on the Tay- 
lor method in Refs. T^and'S and is outlined in Appendix 
IXI Such an algorithm is extremely useful as a check for 
nonstationary problems in which the integration time can 
be very long. The new algorithm agrees with the analyt- 
ical result up to time steps of size, St — .5, in the dimen- 
sionless units of time defined in Eq. |7| In general, the 
time step required to achieve a given accuracy decreases 
as either the frequency or friction increases. 

These results are limited to diffusion in one dimension. 
When studying the motion in the plane defined by the 
rotating magnetic field an average has to be taken over all 
the directions. The integrated friction function, Eq. 1111 
for a particle with the initial velocity inclined with angle 
LOT relative to the magnetic field at the time, t = 0, is 



G{t + T) = 70 a'^t+—sm{uj{t + T)) 
I ^ 
t sin[2w(t + T) 

^2 ^ Aw\ 



(14) 



After some elementary algebra, the integration in ll3l with 
G as in Eq. E| for the case of a constant magnetic field 
leads to the average diffusion coefhcient, 
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70/9 [a? - 1)3/2 



(15) 



of a Brownian particle in a plane. The diffusivity of the 
Brownian particle in a rotating field at various frequen- 
cies has been obtained numerically and is shown in Figure 
13 As can be seen, the diffusivity is a monotonic decreas- 
ing function of the frequency. This result suggests the 
use of the applied field frequency to control the diffusive 
transport of the Brownian particle. 

In the a = 1 limit of this model, there is a divergence in 
the averaged diffusion constant over all the directions at 
constant magnetic field. This limiting behavior is a con- 
sequence of a transition from diffusive to ballistic motion 
at the inclination in which 9 ~ -k. That is to say, that 
it is an artifact of the model in so far as the physical 
system it represents would never take on the value of 
a = 1, and hence would not exhibit an infinite diffusion! 
Nonetheless, the model above serves to demonstrate the 
the accuracy of the numerical and analytical formalism 
when a is far from 1. 



V. REACTION RATES IN A TIME-PERIODIC 
FRICTION 

The assumptions introduced in Sec.^are also applica- 
ble to the description of the reactive interactions between 
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FIG. 3: The normalized average displacement of an ensemble 
of free Brownian particles in the presence of a periodic friction 
is displayed as a function of time. The driving frequencies are 
Lo = 10, 1, 0.1, 0.2 and 0. The result for the fixed case = 0) 
has been calculated analytically using Eq. 1151 The remaining 
results are obtained numerically by averaging over Brownian 
particles with velocities in random inclination relative to the 
magnetic field at the initial time, t = 0. Note that the slopes 
— viz., the diffusion rate — increase with decreasing uj. 



two Brownian particles. Neglecting the hydrodynamic in- 
teraction as before, the dynamics can be described by the 
time-dependent Langevin equation. 



-^V{q) - joHt)q + VXt)C(t) 



(16) 



where q is now a relative mass-weighted coordinate be- 
tween the interacting particles, and V is the potential of 
mean force between the particles. The remaining symbols 
are the same as in the previous section. Phenomenologi- 
cal rate constants — e.g., transition from one metastable 
state of the potential to another or to infinity — cannot 
be calculated analytically when the potential is of a more 
complex form than that of the harmonic oscillator. Di- 
rect numerical evaluation of these rates is usually quite 
time consuming because the time scales involved in the 
problem are widely varying. The reactive flux method 
reduces much of the computation time by initiating the 
trajectories at the barrier.^ It has been used to obtain re- 
active exact thermal escape rates in the stationary limit 
both numerically and exactly, and to obtain approximate 
rates under a variety of limiting approximations. In the 
present case, the problem is nonstationary at short times 
but retains an average stationarity at sufficiently long 
times. The strategy is consequently to generalize the 
rate formula for stationary systems. It must now include 
processes in which stationarity is required only when the 
observables are integrated over a period equal to that of 
the external periodic perturbation. 

In all of the calculations performed here to illustrate 
the approach, the potential has been chosen to have the 
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the form of a symmetric quartic potential, 

V{q) = - 2g2 , (17) 

in which the two minima represent two distinct 
metastable states separated by a dimensionless barrier 
of unit height. (Note that for simphcity, all observables 
in this work are written in dimensionless units relative to 
the choice of this effective barrier and the particle mass.) 
The reactive rate has been calculated for particles with 
inverse temperatures, (3 — 10, or 20. These temperatures 
are low enough to give a well-defined phenomenological 
rate when the reactive flux method is employed in the 
constant friction case, but not so low that trajectories 
are needlessly slow even when one obtains the rate by 
direct methods. 



A. Rate Formula and Numerical Methods 

The standard approach for calculating reaction rates, 
"the reactive flux method," assumes stationarity.24^ In 
establishing its validity, the rate formula needs to be 
checked by comparison with direct methods measuring 
the phenomenological rates between reactants and prod- 
ucts. In this section, a direct approach for obtaining 
rates in the nonstationary cases of interest to this work 
is reviewed and similarly validated. The results of this 
approach are subsequently used to motivate an averaged 
reactive flux formula appropriate for the nonstationary 
case. 

In the direct approach, one simply calculates the rate 
of population transfer from the reactant population ria 
to the the product population n^. The initial popula- 
tion is assumed to be thermally distributed entirely at 
the reactant side. The latter assumption is valid because 
the Boltzmann distribution is the steady state solution 
of the system restricted to the reactant region fApp. IB|| . 
Assuming that a simple first-order master equation de- 
scribes the rate process fApp. IbIi. the population in the 
reactant well can be solved directly as. 



"•a(i) 



«a(io) 



exp 



dt'\'{i) 



to 



(18) 



where the relaxation rate X{t) = k'^ + k~ is the sum of 
the forward (k^) and reverse (fc~) rates, fia is the popu- 
lation in the left well at equilibrium. At equilibrium, for 
a symmetric potential, na(i) = nc{t) = N/2, where N is 
the total population of Brownian particles. In a nonequi- 
librium bath, such as is seen in the model described in 
Sec.m that has oscillatory components with a maximum 
recurrence time, then a phenomenological rate may still 
be obtained by averaging at sufficiently long times com- 
pared to the maximum recurrence time. In particular. 



J to 



\{t')dt' 



■In 



n^{t)-N/2 
N/2 



(19a) 
(19b) 



The second equality was introduced by PoUak and 
FrishmanM as a construction that can lead to long time 
stability thereby ensuring a substantial plateau time4^ 
The instantaneous flux can be found using the differen- 
tial expression^: 



1 



T^a — nc dt 

d 



(na - "c) 



dt 



In (rta - "c) 



(20a) 
(20b) 



The numerical calculation of either of the direct rate 
formulas requires the direct integration of a large number 
of trajectories all initiated in the reactant region. Con- 
sequently, it will only be accurate when the numerical 
integrator is accurate for times that are sufficiently long 
to capture the rate process. This holds at the moderate 
temperatures (near (3V^ = 10) explored in this work for 
which Eas. 1191 and 1201 lead to the same result. The first 
method was used in all of the direct calculations in this 
work because it tends to be more stable. 

The direct methods are time consuming and it is prac- 
tically impossible to apply them at low temperatures. As 
was mentioned at the beginning of this section, the typ- 
ical solution of this problem is the use of the reactive 
flux method. It samples only those states that traverse 
the dividing surface. For stationary systems, the reactive 
flux is- 



(21) 



where the characteristic equation 0[<7(t)] for a trajectory 
is 1 if the particle is in the right well at time t and zero 
otherwise, and the Dirac 5-function ensures that all the 
particles are initially located at the barrier (at a; = 0). 
The angle brackets represent the averaging over the ther- 
mal distribution of the initial conditions. 

One might naively assume that the rate expression 
in Eq. [2] might still hold in the nonstationary case of 
time periodic friction, Eq. |S| The direct and reactive- 
flux rates at different frequencies and different friction 
constants are compared in Table HI The two don't al- 
ways agree and the difference can be as much as an order 
of magnitude. This result should not be surprising be- 
cause of the nonstationarity of the problem. However, 
correlation functions for this system do become station- 
ary when one averages over the period of the external 
perturbation. This suggests that Eq. |^ should be fur- 
ther averaged over the initial time during a period of the 
external field, yielding the average reactive flux rate. 



zvr 



dr 



{S[x{T)]±{T)e[x{t + r)]) 



(22) 
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Rates at 7 = 0.05 
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TABLE I: The integral method of Eq.^|is compared to the 
stationary reactive flux method of Eq. 1211 in calculating the 
activated rate across the double-well potential in a rotating 
field of frequency u>. All the calculations are performed at the 
same bath temperature such that (3V^ = 10, and at three dif- 
ferent values of 7 illustrative of the low, intermediate and high 
friction limits. Here and elsewhere, all values are reported in 
the dimensionless units of Eq. 1161 



(There is a formal proof in Appendix^) A comparison 
between the direct rates and the average reactive flux rate 
is presented in Table ^ The numerics were performed 
at a temperature {(3V^ = 10) that is high enough to 
enable direct calculation of the rate within a few hours 
of CPU time on a current workstation. As can be seen 
from the table, there is very good agreement between the 
methods. Eauationl22lis the central result of this article, 
and represents the fact that the reactive flux method is 
valid for the case of a time- dependent bath when a proper 
averaging is taken over the period of the external field. 
This result is critical for the numerical calculation of rates 
because the direct approaches are cost prohibitive when 
the temperature is much smaller then barrier height. In 
this section the reactive flux method has been generalized 
to include out-of-equilibrium systems in cases in which 
an external force perturbs a bath that is coupled to a 
reactive system. The resulting thermal flux is defined 
only after averaging over the time period of the external 
perturbation. Using the non averaged rate expression 
would lead to undefined rates because the reactive system 
is so far out of equilibrium. A detailed explanation can 
be found in App. ^ 



B. Analytical Approximations 

1. Weak Friction 

For the stationary problem, Melnikov and Meshko\ii^ 
developed a perturbation technique to find the reactive 
fiux at weak to moderate friction limit The expansion 
parameter of the method is the energy loss 6 that a par- 
ticle starting at the barrier experiences while traversing 



the well. Its value is 

<5 = 7/?s , (23) 

where s is the action of a frictionless particle starting 
with zero momentum at the barrier and returning back 
to the top of the barrier after traversing a periodic orbit 
(the instanton), i.e., 

/oo j-qioc) 
p^{t)dt ^ pdq. (24) 

-00 q( — 00) 

The resulting rate is 

k = A:tstT , (25) 

where /ctst — the transition state theory rate 

{u! is the frequency at the bottom of the reactant well and 
is the barrier height) and the depopulation factor T 

is 

1 



1 _ e*(^ +1/4) 



A2 + 1/4 



dX U26) 



The nonstationary analytic rate expression can now 
be obtained by analogy to the formulation of the average 
reactive fiux rate in which the rate is averaged over the 
period of driving term. In particular, the energy loss. 



S{T)=-fl3 ij{t + TYp^{t)dt , 



(27) 



is now obtained as a function of the possible initial con- 
figurations of the driving term which are, in turn, pa- 
rameterized by the time lag r relative to the start of an 
oscillation in the friction of Eq. |S1 Trajectories in one 
dimension can be calculated up to a quadrature directly 
from energy conservation. 



q{to) + 2 f dt ^E^V{q) 

Jto 



(28) 



In the case of the double- well potential defined by Ea.lTTI 
the instanton a,t E = — viz. the periodic orbit on 
the upside-down potential — can be obtained analytically. 
The results for time, 



and momentum. 



^ln( ^+^^ 



(29) 



p{q) = V2q^2 - q^ , (30) 

as a function of q follow readily. By substitution into 
Eq. |57| the energy loss parameter is obtained directly 
with respect to the time lag r relative to the start of an 
oscillation in the friction of Eq. |H1 i.e., 

Sir) = 27/3 [ dqla 



cos 



■^V2qV2 




+ UJT 



(31) 
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TABLE II: The integral method of Eq. 1191 is compared to the average reactive flux method of Eq. 1221 in calculating the 
activated rate across the double-well potential in a rotating field of frequency tJ and various friction constants 7. The potential 
and inverse temperature (/JV* — 10) are the same as in Table |I] At each entry, the integral method result is written above the 
more approximate average reactive flux result. To aid the eye, the latter is also signaled by parentheses. 



The nonstationary rate formula for a time-periodic driv- 
ing friction can thus be written as the product of the TST 
rate and a generaUzed depopulation factor averaged over 
r, 



T[5] 



dTT((5(T)) 



(32) 



in analogy with Eq. 1251 

The validity of the analytical result of Eq. |221 for the 
rate can be checked in the low friction regime in which 
T((5) « (5. Taking the average over a period yields the 
result. 



T[5] 



2^ 



6{T)dT 



3'" 



(33) 



This result is in good agreement with the averaged re- 
active flux rate formula of Eq. as shown in Table IIIII 
at a low friction value (70 — .005), — 10, and var- 
ious frequencies. Even within this weak friction regime, 
as the friction increases, the approximation leading to 
Eg. 1331 will break down. The direct evaluation of Ea. 1321 
corrects this error, and also leads the rate to depend on 
the frequency of the driven friction. 



2. Strong Friction 

The reaction rate in the overdamped regime of the LE 
is well knowni^ The central idea is that the motion in 
phase space is strongly diffusive in this regime. The rate 
is consequently dominated by the dynamics close to the 
barrier. At the vicinity of the barrier top, the potential 
can be approximated by an inverted parabolic potential 
and the LE at the barrier can be written as 



-ujlq - 79 + ^(t) 



(34) 



Rates at 7 = 0.005 



0.1 



10 



MM (Eq.OSJ 8.54x10"^ 8.54x10"^ 8.54x10"^ 
k{t) (Eq.|m 7.5 xlO"^ 7 xl0~^ 7.5 xlO"^ 

TABLE III: The transmission coefficients for the escape rate 
k across a quartic potential at (3V^ — 10 and 7 — .005 ob- 
tained using the average reactive flux method of Eg. 1221 with 
the analytical Melnikov-Meshkov expression 1331 for S. An en- 
semble of 100,000 trajectories has been propagated in each of 
the reactive flux calculations. 



where the fixed friction 7 and stochastic force ^ (t) satisfy 
the regular fluctuation dissipation relation (Eq. ^ and 
iujb is the imaginary frequency at the barrier. It was 
shown that the reaction rate for this case is^L^ 



Ab Wo 
[jjb 27r 



exp 



(35) 



where loq is the frequency of the reactant well, and iAb 
is the imaginary eigenvalue of the homogeneous part of 
Eg 1341 The latter is related to the exponential divergence 
in the trajectories near the barrier. 



q{t) oc e 



(36) 



At strong friction in the nonstationary problem, the 
reaction rate expression is also dominated by the tra- 
jectories in the barrier region. Eguation 1351 can still be 
used for the rates, though now X{t) is the time-dependent 
eigenvalue of the homogeneous part. 







(37) 



of the nonstationary stochastic equation of motion. The 
solution of this equation is not trivial. A possible way 
to solve the problem is found in Ref. |^ It is easier to 
extract the eigenvalue numerically from the exponential 
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FIG. 4: The activated escape rates k of particles in a quartic 
potential and solvated by an anisotropic time-dependent liq- 
uid is obtained as a function of the driving frequency ui using 
two numerical methods described in this work. The numer- 
ical direct rate of Eq. 1191 is shown by dashed lines, and the 
averaged reactive flux rate of Eg. 1221 is shown by solid lines. 
In the former, an ensemble of 250,000 initial conditions were 
used to achieve convergence, and the corresponding numerical 
values are summarized in Table UTTI In the latter, the aver- 
age was performed over an ensemble of 100,000 trajectories, 
yielding the results in a wall-clock time that was an order of 
magnitude faster than that for the direct rate calculations. In 
all cases, the inverse temperature /3V* is 10. 



divergence of trajectories starting near the barrier top, 
g(t)cxe.''*^^(*')'^*' . (38) 

The periodicity of the time dependent coefRcient in 
Eq. 123 leads also to a periodicity in Ab(i). If Ab is the 
time average of the time-dependent eigenvalue of Eq. 1371 
over a period, then for t much larger than the period, 
Eq. |2Hlis analogous Eq. [3^1 with Ab in the exponent. Re- 
placement of Ab by Ab in the rate expression (Eq. I35() 
provides good agreement with the averaged reactive flux 
rates as shown in the high friction columns of Table IIVI 

3. Weak to High Friction 

The results of the two previous subsections have mo- 
tivated the redefinition of the components of the rate 
formula in the low and high friction limits of the nonsta- 
tionary time-periodic problem. Retaining these assign- 
ments in the stationary turnover rate formulaSS suggests 
the nonstationary turnover rate. 



k = 



Ab (my 

tJb 27r T[25] 



■ exp 



(39) 



FIG. 5: The activated escape rates k of particles in a quar- 
tic potential and solvated by an anisotropic time-dependent 
liquid is obtained as a function of the driving frequency u) 
comparing the reactive flux approach to an analytical result. 
(The latter is expected to be accurate here — and not in Fig. |1| 
or Table El — because the inverse temperature has been in- 
creased to 20.) The solid line corresponds to the numerical 
result calculated using the averaged reactive flux rate method 
of Eg. 1221 and the dashed line corresponds to the turnover for- 
mula in Eq. 1391 The numerical calculations were performed 
by averaging over an ensemble of 250,000 trajectories. 
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The prefactors from the nonstationary turnover rate are 
compared to those from the averaged reactive flux rate at 



10.0 (.155) (.727) (.82) (.685) (.31) (.21) 
.148 .720 .821 .683 .3 .219 

TABLE IV: The average reactive flux method of Eq. is 
compared to the analytic approximation of Eq. 1391 in calcu- 
lating the activated rate across the quartic potential in a ro- 
tating field of frequency u and various friction constants 7. 
The inverse temperature [PV^ = 20) is higher in contrast to 
the previous tables. At each entry, the more approximate av- 
erage reactive flux result is written above the analytic result. 
To aid the eye, the former is also signaled by parentheses. 



the inverse temperature /3 = 20 in Table Hvl and in Fig.El 
As can be seen, there is a very good agreement between 
the numerical and analytic results at the very weak and 
strong friction limits. Therein the results diff'er by no 
more than 5% throughout the frequency range; an er- 
ror margin smaller than the error bars in the numerical 
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calculations. At moderate friction and low frequencies, 
however, the differences — on the order of 10% — cannot 
be explained by error in the numerical calculations alone, 
and may be significant. Corrections or improvements in 
the approximations leading to the connection formula of 
Eq.|2niare also of interest, but not pursued further in this 
work. Recall that the turnover escape rate expression for 
the LE with constant friction was obtained through the 
solution of the equivalent Hamiltonian formalismi^ A 
similar approach for the solution of a the Hamiltonian 
equivalent of the stochastic time-dependent bath prob- 
lem may lead to a fruitful solution. However, even in the 
constant friction case, the turnover formula can give rise 
to small systematic error. With these reservations, the 
approximate rate formula can be used to obtain time- 
dependent escape rates. It is clear from the results that 
there is a frequency effect on the reaction rates. For the 
specific example studied here, the effect can modify the 
reaction rate 



VI. DISCUSSION AND CONCLUSIONS 

In this work, several techniques for obtaining the dy- 
namics of interacting Brownian particles that are coupled 
to a time dependent thermal bath have been discussed. 
Two models, one of dynamics in lyotropic liquids and 
one for dynamics in pure nematic liquid under a periodic 
external field has been brought as examples of such sys- 
tems. The models include a new mechanism for stochas- 
tic dynamics in which an external force is used to drive 
the thermal bath. There is no net injection of energy 
to the Brownian particles in the bath due to the driv- 
ing force, hence they keep their equilibrium properties. 
Yet observables such as reaction and diffusion rates are 
modified. The existence of a steady state that retains the 
equilibrium enables one to express out-of-equilibrium ob- 
servables with respect to averaging over the equilibrium. 
This is the Onsager regression hypothesis ( Appendix 0. 
We used this to extend known methods for calculating 
the reaction rates in the constant friction to nonstation- 
ary baths. Extensive computation effort was used to il- 
lustrate the diffusive and reactive rates for an effective 
Brownian particle in the naive anisotropic liquid bath 
model with rotating magnetic field. However, the nu- 
merical and analytical tools that have been modified and 
developed in this work are appropriate for any model 
with time-dependent friction. 

The construction introduces new control parameters 
into the problem; namely, the external force amplitude 
and frequency. We concentrated on the latter and exhib- 
ited the frequency dependence of diffusion and reaction 
rates in the naive model. This dependence is not lin- 
ear and changes dramatically with the friction strength. 
The enhancements in the reaction rate and the diffusion 
coefficient are not the same, i.e., the maximum in the 
diffusion rate as a function of the external frequency is 
not the same as the maximum in the rate. This nonlin- 



ear behavior could be used to enhance reaction diffusion 
processes, such as cluster nucleation, by up to a few order 
of magnitude. 

In the extension of the naive model to more realistic 
nematic liquids, the cooperative effects of the liquid can- 
not be omitted. There are phenomenological difficulties 
in defining friction and the fluctuation dissipation rela- 
tion in liquid crystals. To the best of our knowledge such 
a theory is still not fully developed. The development 
of such a theory based on microscopic assumptions is 
an extremely challenging problem. Boundary effects and 
elastic forces will create dynamical micro-domains char- 
acterized by differing uniform directors in a real nematic 
under rotating magnetic field. The theory for dynamics 
in nematics will have to deal also with the spatial inho- 
mogeneities. These are among the challenges to future 
work in trying to better understand the diffusive dynam- 
ics in lyotropic liquids. 
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APPENDIX A: FOURTH-ORDER INTEGRATOR 
FOR THE LE WITH PERIODIC FRICTION 

A high-order integrator was developed for the regu- 
lar LE or CLE in Ref. 42. A modified algorithm for 
time- and space-dependent friction was developed for 
the explicit GLE with exponential memory kernel in 
the frictioni^. This appendix introduces the numerical 
scheme necessary for solving a time dependent stochastic 
equation equation of motion of the form of Eq. [71 

A finite difference scheme is used to propagate the solu- 
tion over a small time step. At each iteration, the prop- 
agator is expanded to fourth order with respect to the 
time step using a strong Taylor scheme The resulting 
integrator can be decomposed into two uncoupled terms: 

q{t + h) = qdct{p,q,t) + qran{p,q,t) (Ala) 

p{t + h) = pdct{p,q,t) +pran{p,q,t) . (Alb) 

The deterministic terms that are collected within gdot 
and pdct are those that remain in a fourth-order Taylor 
expansion of the deterministic equation of motion after 
removing any term that includes a stochastic variable. 
The deterministic propagator can be calculated numer- 
ically with any fourth-order deterministic scheme; the 
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fourth order Runge Kutta method was chosen for this 
work.'*® 

The random propagator (giving rise to the stochas- 
tic contribution, gj-an and Pran) includes all terms up to 
fourth-order that include stochastic variables. For the 
present case, the random integrator for the space and 
momenta leads to the contributions. 



+ + -i!h^ 



-V"^j + -f^^/ + 2i; + 7-/^^i; Zij (A2) 



■0 + -0 + ^-iph^ 



-7-0^ - -0 - (Vj + 3^ip'^ip)h 
- (-^^■(')+7(3^^' + ^^'^))/i' 



v"ip + j'^tp^ + V' + -iji^^-^ 

+ {p°v"'ip + 57V^^ + 47(21/'^/'^ + V'^V') 



,(3) 



2jV -j^^p^ +poV V + i^-^j 
- 97V*^ - 47(2V'^/'2 -f V'^VS) Zij , (A3) 



where ip{t) is the coefficient in Eq. |S1 dictating the time- 
dependence in the friction, and the Gaussian random 
variables {Zi} are correlated according to the moments 
specified in Ref. li^. The symbol, with the paren- 

thesized superscript, denotes the third-order derivative 
of tp with respect to time. There is frequent repetition of 
terms differing by no more than a ratio of constant co- 
efficients, and this allows the integrator to be computed 
very economically despite the seemingly large number of 
terms contained above. In fact, the most time-consuming 
part of the scheme is the calculation of the random num- 
bers. Although the fourth-order integrator has been ex- 
panded in terms of the stochastic variables, the neglected 
higher-order terms have coefficients that depend on the 
friction and frequencies to high order. Consequently, the 
algorithm loses its efficiency at high friction or in the 
high frequency domain. For a given problem, a compar- 
ison of the fourth-order algorithm to lower-order algo- 
rithms such as the velocity Verlet algorithm*^ is advis- 
able to ensure that the requisite accuracy is achieved. 
Though not presented explicitly here, such comparisons 
have been performed with favorable agreement for the 
models of this work. The fourth-order integrator is ac- 
curate, and equally importantly, provides a substantial 
savings in computational time. 



APPENDIX B: REACTIVE FLUX METHOD FOR 
SYSTEMS IN NONSTATIONARY 
ENVIRONMENT 



The reactive flux formalism is an efficient numerical 
tool for calculating the escape rates of Brownian par- 
ticle from a metastable well at low temperatures. The 
derivation of the reactive flux method is well formulated^ 
In this appendix, the derivation is recapitulated in oder 
to emphasize the new considerations that emerge be- 
cause of nonstationarity. The corner stone for any rate 
calculation in a nonequilibrium system is "the Onsager 
regression hypothesis."— This hypothesis asserts that 
"the relaxation of macroscopic non-equilibrium distur- 
bances is governed by the same laws as the regression of 
spontaneous microscopic fluctuations in an equilibrium 
system."— The two basic assumptions of the regression 
hypothesis are the existence of an equilibrium and a con- 
siderable (large) separation of time scales between the 
motion of the subsystem and that of the overall relax- 
ation of the system. 

The proof of the existence of an equilibrium in a sys- 
tem is more readily obtained through an analysis of the 
probability distribution W{q,p;t) of the Brownian parti- 
cles rather than the explicit trajectories calculated using 
the LE. The equation of motion for the distribution is 
the Fokker-Planck equation, i 



d_ 

dt 



W{q,p;t) 



dq dp 
dp^ (i 



2 y{t)\W{q,p-t) 



(Bl) 



which corresponds to the nonstationary LE in Eq. [7| A 
direct substitution of the Boltzmann distribution con- 
firms that it is indeed a solution of this equation. From 
this one can deduce that the Boltzmann distribution 
is the equilibrium distribution for systems with time- 
dependent friction. It might seems strange that equi- 
librium doesn't change when the system is driven by a 
time-dependent forcing friction. The key point is that the 
force is coupled only to the bath. The bath is taken to 
be infinite dimensional and absorbed the energy extract 
by the force. This, of course, does not mean that the sys- 
tem dynamics is not modified. On the contrary, as has 
been shown in the article, the diffusivity and reaction 
rates are modified by the periodicity of the externally 
driven friction. The time scales relevant to the second 
assumption are the period rt of the external force, the 
transient time Tg of the system, the typical escape time 
To of a thermal particle and the observation time r. By 
construction of the naive model and the choice of its pa- 
rameterization, these times follow the simple inequality: 
(rf or Ts) <C r <C To. As such, the system necessar- 
ily satisfies both of the assumptions needed to apply the 
Onsager hypothesis. 

The rate constants may be obtained from a first order 
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master equation representing the population dynamics, 

Mt) = -k+{t)nS) + k-{t)n^{t) (B2a) 
n,{t) = k+{t)n,(t) - k-{t)n,{t) , (B2b) 

where {ric) is the population of the left (right) well of 
the quartic potential in Eq. 1171 The nonstandard feature 
is that the forward (backward) coefhcient rate k^{t) is 
time dependent. As described in the text, they do have 
the same periodicity as the time-dependent friction The 
master equation in Ea. lB2l can be integrated in the usual 
manner to obtain the result. 



Ana(t) 
Ana(to) 



-/,* \{t')dt' 
e 'o 



(B3) 



where Ana(t) is the fluctuation of the momentary popu- 
lation of the left well relative to the equilibrium popula- 
tion of the well, and A = (fc+ + fc^) as in the text. The 
Onsager regression hypothesis enables the connection to 
the equilibrium averaged expression, 



m<i{h)?) 

Taking the time derivative of both sides leads to 



(B4) 



{5e[q{t^)]6e[q{t)] 



^A(t)e"^*o^(*')''*' (B5a) 
^A(t) . (B5b) 



The last equality is a result of the large time scale separa- 
tion. So far there is no real difference from the standard 



derivation. In the usual derivation, the assumption of 
stationarity would now permit the modification of the 
numerator into the form of a flux correlation function. 
This cannot be performed in the present time-dependent 
case. However, stationarity can be regained in this sys- 
tem by averaging over the period of the time-dependent 
friction. In practice, this can be achieved by initiating 
the subsystem at various time shifts r relative to some 
arbitrary time origin in the time-dependent friction. An 
integration over the period of IB5I leads to the averaged 
form of the reactive flux. 



dr 



27r 



(B6a) 







2^ JO 

-A 



mq{T)Y) 
dT\(T) (B6b) 
(B6c) 



The bar in the definition of A indicates the average over 
time period of the time-dependent friction. Equation IB6I 
is the averaged reactive flux rate for the nonstationary 
system with periodic friction. The above algebra also 
justifles the extension of the Melnikov-Meshkov theory to 
the time dependent friction case averaged over a period 
of the external perturbation described in Sec. IV Bl 
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